Influence of self-gravity on the runaway instability of black hole-torus systems 
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Results from the first fully general relativistic numerical simulations in axisymmetry of a system 
formed by a black hole surrounded by a self-gravitating torus in equilibrium are presented, aiming 
to assess the influence of the torus self-gravity on the onset of the runaway instability. We consider 
several models with varying torus-to-black hole mass ratio and angular momentum distribution 
orbiting in equilibrium around a non-rotating black hole. The tori are perturbed to induce the mass 
transfer towards the black hole. Our numerical simulations show that all models exhibit a persistent 
phase of axisymmetric oscillations around their equilibria for several dynamical timescales without 
the appearance of the runaway instability, indicating that the self-gravity of the torus does not play 
a critical role favoring the onset of the instability, at least during the first few dynamical timescales. 
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Introduction. Self-gravitating tori orbiting black holes 
(BHs) may form after the merger of a binary system 
formed by a BH and a neutron star (NS) or the system 
formed by two NS (see e.g.[l[ and references therein). In 
addition, they may also be the result of the gravitational 
collapse of the rotating core of massive stars 0, Q ■ State- 
of-the-art numerical simulations have started to provide 
quantitative estimates of the viability of such systems to 
form [l|, As such BH-torus systems are thought 

to be the central engine for gamma-ray bursts (GRBs) 
12, 13]. understanding its formation, dynamics and sta- 
bility properties is of high relevance. 

In particular, the so-called runaway instability, first 
found by Abramowicz, Calvani and Nobili [T3], is an ax- 
isymmetric instability that could destroy the torus on 
dynamical timescales. In a marginally stable torus, the 
radial pressure gradient may drive the transfer of mass 
towards the BH through the cusp-like inner edge of the 
torus. Due to the accretion of mass and angular mo- 
mentum, both the mass of the BH and its spin increase, 
and the gravitational field changes, leading to two possi- 
ble evolutions: (i) if the cusp moves inwards towards the 
BH, the mass transfer slows down and the system is sta- 
ble, or (ii) if the cusp moves deeper into the torus, mass 
accretion will increase, and the accretion process will be 
runaway unstable. 

The numerical study of the runaway instability has 
so far been investigated under different approximations 
(see e.g. [IH, HH). Abramowicz et al. [lj], assuming a 
pseudo-Newtonian potential for the BH, constant angular 
momentum distribution in the torus and an approximate 
treatment of the disk's self-gravity, found that the insta- 
bility occurs for a wide range of initial models. More 
detailed studies based on stationary models, either as- 
suming a pseudo-Newtonian potential for the BH |l3] or 
being fully relativistic calculations [l8|], indicated that 
the self-gravity of the disk favors the instability, by ar- 
guing that, as a result of the accretion process, the cusp 
would move closer to the center of the torus than in non 



self-gravitating disks. However, there are additional pa- 
rameters which have a stabilizing effect: (i) the rotation 
of the BH [19( , and (ii) the most important one, the ra- 
dial distribution of specific angular momentum, increas- 
ing with the radial distance [201 ] . 

The first time-dependent, general relativistic hydrody- 
namical (GRH) axisymmetric simulations of the run away 
instability were performed by Font and Daigne (l5l Il6j . 
The BH evolution was assumed to follow a sequence of 
stationary BH spacetimes of increasing mass and angu- 
lar momentum, controlled by the mass and angular mo- 
mentum transferred from the torus, whose self-gravity 
was neglected. The first work [lj| , which focused on tori 
with constant distribution of specific angular momentum 
/ = —u v /u t , with u v and u t being the corresponding 
components of the 4- velocity u^, showed that the sys- 
tem is runaway unstable on a dynamical timescale. On 
the other hand, the second work [l|| showed that thick 
disks with non-constant specific angular momentum dis- 
tributions, increasing outwards with the radial distance 
according to a power law I = Kr a , are stable for very 
small values of the angular momentum slope a (much 
smaller than the Keplerian value a = 0.5), confirming 
the prediction of stationary studies. 

Despite the progress that has been made the existing 
works are still not conclusive, mainly due to the absence 
of important physics in the modeling. The complexity of 
handling the presence of a spacetime singularity in addi- 
tion to the hydrodynamics and the self-gravity of the ac- 
cretion torus, make full GRH simulations of such systems 
very challenging. The simulations presented in this Letter 
accomplish the goal of assembling this important physics 
to provide a conclusive answer about the likelihood of the 
instability on the first few dynamical timescales, which 
cause the most concern in the context of models of short 
GRB. Although we do not exclude the onset of the insta- 
bility on much longer timescales, it is only meaningful in 
connection with other relevant instabilities, specially the 
magnetorotational instability. 
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TABLE I: Main properties of the equilibrium models studied in units of c = G — Mq = 1 (unless shown otherwise). From left 
to right the columns show: the type of specific angular momentum distribution, the torus-to-BH mass ratio, the position of 
the maximum density point r max , the position of the inner and outer radii of the torus r- m and r out , the maximum rest-mass 
density and the orbital period at the center of the torus t or b- 

Model j - law M t /Mbh ^max ?"in Tout Pmax (ff/cm 3 ) t orb 



Ml const 0.1 7.17 4.92 10.17 3.189 x 10 14 147.81 

M2 const 1.0 8.87 4.02 19.97 2.202 x 10 14 199.54 

M3 non-const 0.1 10.47 4.92 19.97 3.902 x 10 13 245.37 

M4 non-const 0.5 10.02 4.07 19.97 1.538 x 10 14 229.91 



Numerical setup. The numerical simulations have been 
carried out with the nada code (see [2l| for details and 
tests) which solves the Einstein equations (using the 
BSSN approach, the "cartoon" method, and the mov- 
ing puncture approach) coupled to the GRH equations 
(solved with a high-resolution shock-capturing scheme 
based on 3rd-order PPM reconstruction and the HLLE 
approximate Riemann solver) . In addition, we use a stan- 
dard T-law equation of state (EOS), for which the pres- 
sure is expressed as a function of the rest-mass density 
and specific internal energy as P = pe{T — 1), where T 
is the adiabatic exponent. For the vacuum region out- 
side the torus we use a dynamically unimportant artifi- 
cial atmosphere with rest-mass density pthr ~ 10~ 8 /9 max - 
The evolution equations are integrated by the method of 
lines, for which we use a 4th-order Runge-Kutta scheme. 
For the simulations reported here we use an equidis- 
tantly spaced (x, z) grid with a grid spacing Ax = Az = 
0.05M B h and N x x N z = 600 x 600 points to cover a 
computational domain, < x < L and < z < L, with 
L = 30Mbh- Unless otherwise stated we use units in 
which c = G = M Q = 1. 

Initial data. Compact equilibrium configurations for a 
BH-torus system are obtained in the moving puncture 
framework (we refer to [22j for details). We adopt a 
r = 4/3 polytropic EOS to mimic a degenerate relativis- 
tic electron gas, and construct initial configurations using 
either constant or non-constant specific angular momen- 
tum distributions, defined as j = hu v (h being the spe- 
cific enthalpy). A list of the models along with their main 
features is given in Table [U We consider four different tori 
around a non-rotating BH (of mass Mbh = Mq)- Follow- 
ing [H, [H, HI , the EOS polytropic constant n is chosen 
such that the torus-to-BH mass ratio, Mt/MeH, is 0.1, 
0.5 or 1, depending on the model. We note that existing 
simulations of NS/NS and BH/NS mergers [J EMM if 
yield M t ~ 0.01 — 0.2 M Q . As these configurations are 
not overflowing the cusp, we introduce an initial pertur- 
bation to induce a small mass transfer through the inner 
edge of the tori (as in In our simulations we simply 

perturb the v x component of the 3- velocity of the torus 
as v x w —77, which otherwise would initially be zero. 
For each model the numerical simulations are stopped at 
t ~ 2000, which corresponds to ~ 10 ms (between 8 to 
10 orbital periods depending on the model), since at late 



times the growth of the Hamiltonian constraint violation 
would lead to a loss of accuracy for the spacetime evo- 
lution. Nevertheless, the time scale we are considering 
in our simulations would allow to identify the runaway 
instability, if present, as it could even take place within 
one orbital period for the more massive models M2 and 
M4 (see [2 El]). 

Results. The left panel in Figure [1] displays the time 
evolution of the total rest-mass and central rest-mass 
density, each of them normalized to its initial value, for 
the evolution of models Ml (solid and dashed lines) and 
M2 (dotted line). Ml is a model with j-constant and 
with an initial rest-mass of M t = O.IMbh, thus repre- 
senting a model with torus-to-BH mass ratio in agree- 
ment with results obtained by simulations of NS/NS or 
BH/NS mergers 0, i, 0, S El • For this initial model, 
we have considered two different initial perturbations, 
r/ = 0.01 (solid line) and r] = 0.025 (dashed line) to eval- 
uate its influence on the overall dynamics. As expected, 
the initial perturbation triggers a phase of axisymmetric 
oscillations of the torus around its equilibrium which are 
present throughout the simulation. Such oscillations in- 
duce a small outflow of matter through the cusp towards 
the BH. This, however, does not reduce significantly the 
total rest-mass of the torus, plotted in the upper panel. 
At the end of the simulation (t ~ 10 t rb) the rest-mass of 
the torus Ml is conserved up to about 1%. Therefore, the 
BH mass and spin do not increase with time significantly, 
and the torus shows no sign of the runaway instability. 
Further information about the process of accretion is ob- 
tained from the right panel in Figure [TJ which shows the 
time evolution of the mass accretion rate. We notice that 
the mass accretion rate for model Ml is larger the larger 
the initial perturbation is. For rj = 0.025, there is an 
initial stage of very small mass transfer through the in- 
ner edge of the torus which lasts for about half an orbit 
(t ~ 80). This is followed by a stage in which the oscil- 
latory behavior of the rest-mass accretion rate, signature 
of the induced oscillations, is obvious throughout the nu- 
merical evolution. Interestingly, during the oscillation 
phase, the mass flux does not increase in amplitude with 
time, as one would expect prior to the onset of the run- 
away instability (see [15|, U^). Instead, it reaches a max- 
imum of about M ~ 0.1 Mq/s after the second orbital 
period. Later it decreases and exhibits a series of oscilla- 
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FIG. 1: Left panel: Time evolution of the total rest-mass (top) and central rest-mass density, each of them normalized to its 
initial value, for the evolution of models Ml and M2. Right panel: Mass accretion rate evolution for models Ml and M2. 



tions around a lower value, never showing any signature 
of exponential growth. The mass accretion rates of the 
perturbed tori are in good agreement with the maximum 
expected accretion rates for hyper-accreting disks associ- 
ated with the central engine of gamma-ray bursts, which 
could be as high as Af ~ 0.01 — 1 M Q /s, depending on the 
formation mechanism [25] ■ In the case of a disk formed 
by the merger of a NS with another compact object (ei- 
ther a NS or a BH), most of the material in the accretion 
disk would be accreted on a time scale comparable to the 
viscous time scale t v i sc ~ 0.1 s. In the collapsar scenario, 
such high accretion rates could be sustained for ~ 10 s 
as the disk is fed by the in-falling stellar material. 

We consider next a significantly more massive torus, 
model M2, while keeping the same rotation law. The to- 
tal rest-mass of this torus is M t = I.OMbh- We use a 
small initial perturbation, r\ — 0.01, because due to the 
larger coordinate dimension of the torus, a larger per- 
turbation would cause its outer parts to move outside 
the computational domain. Despite being more massive 
the overall dynamics of M2 is very similar to the one 
discussed above for the less massive torus, Ml. As ex- 
pected, the time evolution of the central rest-mass den- 
sity, displayed in the left panel of Figure Q] with a dot- 
ted line, shows a series of axisymmetric oscillations dur- 
ing the entire length of the simulation. As we have ob- 
served for the low mass torus, the amplitude in the evo- 
lution of the mass flux does not increase with time and 
does not lead to the onset of the instability. We note 
that these results do not actually contradict results ob- 
tained for non-self-gravitating tori with constant distri- 
bution of angular momentum (la ]. Notice that despite 
the difference in the definition of the specific angular mo- 
mentum, the j-constant condition in our models leads to 
Z-constant up to a difference of the order of 10~ 5 . In 




FIG. 2: Mass accretion rate for models M3 and M4. 



addition, models of [l5[ and the ones considered here 
satisfy the condition j'torus > iisco throughout the evo- 
lution, where ISCO stands for innermost stable circular 
orbit. However, those previous studies considered ini- 
tial models which were overflowing the cusp in order to 
induce a large stationary accretion rate, which varied be- 
tween M ~ 0.1 — 34.0 Mq/s. Unstable runaway behavior 
in the case of mass accretion rates of M ~ 0.1 M /s 
were found on a timescale of 100 i or b, and only on a 
dynamical timescale for the largest values of the mass 
flux. Similar results were found by 24] introducing an 
initial perturbation on the equilibrium tori. Since in our 
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case M ~ 10~ 3 — 10~ 4 M Q /s throughout the simula- 
tions, the exponential growth of the mass accretion would 
not manifest itself even on such timescales. This is in 
agreement with the recent work of [l| where a systematic 
study of the BH-torus systems produced by the merger of 
unequal-mass NS binaries was presented, concluding that 
self-gravitating tori with mass accretion rates as high as 
M ~ 2.0 M Q /s were stable on the dynamical timescales 
investigated. 

Motivated by the influence of different rotation laws 
on the onset of the instability [l5| we have also carried 
out numerical evolutions of two j-non-constant models, 
M3 and M4, with two different torus-to-BH mass ratios, 
0.1 and 0.5 respectively. Despite the difference in the 
rotation law with respect to models Ml and M2, the dy- 
namics is very similar, as inferred from the evolution of 
the mass accretion rate displayed in Figure [2] No expo- 
nential growth of the mass flux is found and, therefore, 
no runawa y in stability is present. 

Studies |i~7l . Il8j mainly based on stationary models, 
either assuming a pseudo-Newtonian potential for the 
BH and a Newtonian potential for the self-gravity of 
the torus, or being fully relativistic calculations indicated 
that the self-gravity of the disk favors the instability. The 
simulations presented here overcome the limitations of 
preceding works. Our results indicate that in the general 
case in which mass accretion rates are consistent with 
those expected for hyper-accreting disks, and in which 
the angular momentum distribution increases with the 
radial distance, the effect of self-gravity is not sufficient 
to lead to unstable accretion. 

Despite the simplifying assumptions of our models 
(magnetic fields or detailed microphysics are not consid- 
ered), our results are consistent with the expected dy- 
namics and mass accretion rates for hyper-accreting disks 
and do not challenge the time scale required for produc- 
ing a GRB. Moreover, although magnetic fields have not 
been considered in our simulations, it is likely that these 
do not play a critical role in the context of the runaway 
instability. Despite the fact that numerical simulations 
of magnetized NS mergers are still scarce and the re- 



sults are inconclusive, simulations by [28j indicate that 
the effect of the magnetic fields during the merger and 
the initial phase of the BH-torus lifetime is not dramatic 
(as reflected by the similarities between the gravitational 
waves in the magnetized and unmagnetized cases), rather 
these are expected to become more important during the 
secular evolution of the accretion torus. 

Conclusions. We have presented results from the first 
fully general relativistic numerical simulations in axisym- 
metry of a system formed by a BH surrounded by a 
marginally-stable self-gravitating torus aiming to assess 
the influence of the torus self-gravity on the onset of the 
runaway instability. Several models with different torus- 
to-BH mass ratio and angular momentum distributions 
have been considered. The tori have been perturbed to 
induce mass transfer towards the BH. Our numerical sim- 
ulations show that all models exhibit a persistent phase of 
axisymmetric oscillations around their equilibria for sev- 
eral dynamical timescales without the appearance of the 
runaway instability. Thus, the self-gravity of the torus 
does not play a critical role favoring the onset of the 
instability. Clearly, to investigate additional m = 1 non- 
axisymmetric features that may play a role on the dy- 
namics of the system 3D simulations are required. Nev- 
ertheless, the robustness of our results in axisymmetry 
on the influence of self-gravity on the runaway instabil- 
ity are confirmed by simulations we have also performed 
with the independent 2D code of [26J . These simulations 
and the investigation of non-axisymmetric instabilities of 
self-gravitating tori around BHs will be presented else- 
where. 
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